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1.  Part  I:  Air  Gun  Modeling 


1.1  Introduction 


In  an  air  gun  test,  it  is  critical  to  identify  the  characteristics  of  the  test  environment  such  as  the 
dynamic/physical  properties  of  momentum  exchange  mass  (MEM)  and  mitigator  combinations  in 
order  to  obtain  a  target  forcing  function.  Figure  la  depicts  a  schematic  of  the  air  gun  test  setup. 
The  test  setup  consists  of  an  air  gun,  a  pressure  tube,  a  catch  tube,  and  a  MEM.  The  test  item, 
which  is  placed  inside  the  pressure  tube,  is  launched  to  impact  an  aluminum  honeycomb  mitigator. 
The  test  item  contains  an  on-board  recorder  (OBR)  that  records  responses  during  launch  simula¬ 
tion.  The  OBR  case,  which  is  hereafter  referred  to  as  the  “OBR,”  is  shown  in  figure  lb  with 
several  accelerometers  and  strain  gauges.  The  model  developed  in  this  study  will  be  used  to 
determine  the  response  behavior  of  the  test  item  mounted  on  a  given  projectile  during  a  launch 
simulation  air  gun  test.  This  methodology  requires  the  development  of  a  predictive  model  of 
responses  of  the  test  article.  An  analytical  model  is  herein  developed  to  simulate  an  air  gun  launch 
environment  in  which  a  test  object  mounted  on  a  projectile  is  launched  through  the  air  gun  and 
decelerated  by  the  crushing  of  an  aluminum  honeycomb  mitigator  that  impacts  the  MEM  before 
being  stopped  at  the  retrieving  end. 


CATCH  TUBE 


PRESSURE  TUBE 


GUN  EXIT 


MEM 


1 1  1 

LJ 

MITIGATOR^ 

TEST  Item 


Energy 

Projectile 


a.  schematic  of  air  gun  test  setup. 


b.  The  OBR  with  test  item. 


Figure  1.  Schematic  of  the  air  gun. 
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The  objective  of  this  effort  is  to  develop  an  overall  analytical  model  to  predict  the  transient 
response  of  a  generic  artillery  component  subjected  to  launch  simulation.  The  simulation  is 
achieved  via  impact  mitigation  techniques  whereby  the  component  is  propelled  to  a  target  inside 
a  4-inch  carrier  equipped  to  measure  component  response  during  the  simulation.  A  detailed 
description  of  the  developed  formulation  is  presented  in  the  next  section. 

1.2  Model  Formulation 

1.2.1  Assumptions 

To  model  the  response  of  the  OBR-mitigator-MEM  system  during  impact,  an  analytical  model  is 
herein  developed.  The  model  provides  a  closed  form  solution  to  the  dynamic  system  driven  by 
the  central  difference  equations  and  is  based  on  the  following  assumptions: 

The  OBR  and  MEM  can  be  accurately  represented  by  rigid  bodies  since  they  do  not  dissipate 
energy  through  plastic  deformation.  Only  the  lowest  natural  frequency  of  the  OBR  is  of  interest 
when  we  are  predicting  the  acceleration  output;  therefore,  the  OBR  is  represented  by  a  2-degree- 
of-freedom  spring-mass  system.  It  is  assumed  that  this  frequency  is  readily  available  (if  not,  it 
can  be  easily  detennined  through  simple  analysis). 

The  mass  of  the  mitigator  is  neglected.  The  mitigator  is  entirely  represented  by  the  force  it  exerts 
on  the  OBR  and  the  MEM.  This  force  is  detennined  from  the  geometric  and  material  properties 
of  the  mitigator  and  from  balance  and  conservation  considerations  during  the  analysis  stage. 

1.2.2  Rigid  Body  Dynamics  System 

The  impact  event  is  represented  by  a  rigid  body  mass-spring  system  with  3  degrees  of  freedom 
acting  along  the  horizontal  axis.  The  system  consists  of  three  masses  and  one  spring,  as 
illustrated  in  figure  2.  The  OBR  has  initial  velocity  v0  and  is  represented  by  two  of  the  masses, 
m i  and  mo,  and  a  connecting  spring,  k\.  The  mitigator  is  represented  by  the  force  F,  and  the 
MEM  is  represented  by  its  mass  m3.  Note  that  the  whole  system  is  divided  into  two  separate 
systems  coupled  by  the  force  F  representing  the  mitigator. 

This  assembly  is  a  very  simple  representation  of  a  complex  deforming  system.  If  the  first  OBR 
frequency  mode  is  to  be  neglected,  we  can  further  simplify  it  by  combining  masses  mi  and  m2 
into  a  single  body  of  mass  m0BR  =  mi  +  m2  and  dropping  spring  k\.  This  will  result  in  a  2-degree- 
of-freedom  system.  However,  the  essential  part  of  this  system  formulation  is  the  representation 
of  the  mitigator  with  the  force  F.  To  accurately  represent  the  behavior  of  the  mitigator,  an 
elaborate  and  efficient  formulation  needs  to  be  developed.  Such  a  formulation  is  described  in 
detail  as  follows.  This  simple  dynamic  system  representation  is  accurate,  efficient,  and  provides 
adequate  results. 
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Figure  2.  Dynamic  system  representation. 


First,  let  us  describe  a  time-stepping  scheme  to  solve  the  dynamic  system.  The  system’s  motion 
can  be  represented  by  the  system  of  equations: 

[m}{u}  +  [k]{u  }  —  {f }  (1) 

in  which 


[m]= 

M- 


m]  0 
0  m. 


is  the  mass  matrix  of  the  left-hand  side  system. 


k  j  -  k , 

k  |  kx 


is  the  stiffness  matrix  of  the  left-hand  side  system, 


u 


is  the  displacement  vector  with  ui  being  the  displacement  of  mass  m, 


where  i  =  1,2;  {u}  and  {ii}  are  the  mass  velocities  and  accelerations,  respectively, 
{f }  =  |  ^  1  is  the  external  force  vector. 


To  solve  this  system,  we  define  a  central  difference  time-stepping  scheme.  At  the  initial  moment 
(t  =  0),  we  have 
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Equations  4  and  3  are  sufficient  to  define  a  forward  time-stepping  scheme,  provided  that  a 
sufficiently  small  time  step,  At,  is  used.  To  start  the  time-stepping  scheme,  {u}_i  needs  to  be 
available,  which  can  be  detennined  from  the  second  equation  of  equation  2: 

{u}_!  =  {u}0  -  At{u}0  +  {ii}0 

For  the  right-hand  side  system  representing  the  MEM  we  have 
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m3«3  =  F 

fffT  ((«3  Li  -  2(m3  )n  +  (“3  L,  )  =  K 
At2  F 

(M3  L  =  - -  +  2(«3  )„  -  («3  Ll 

m3 

and  the  velocity  and  acceleration  of  the  MEM  can  be  calculated  from  equation  3. 

1.3  System  Parameters  Definition 

1.3.1  OBR  and  MEM 

The  simplest  part  of  the  system  (figure  2)  is  the  MEM,  which  is  represented  by  its  mass  only. 

The  OBR  is  represented  by  two  masses  and  a  connecting  spring.  The  sum  of  the  two  masses  is 
equal  to  the  total  mass  of  the  actual  OBR.  In  the  determination  of  the  way  in  which  the  total 
mass  is  distributed  between  m\  and  m2,  different  values  were  used  for  the  investigated  system. 

It  was  determined  that  a  ratio  of  m\\m2  anywhere  from  1 :2  to  2: 1  gives  reasonable  results  for  the 
second  mode  or  lowest  fundamental  frequency  of  the  OBR  (the  first  mode  is  a  rigid  body  mode). 
The  spring  stiffness  (k\ )  is  determined  so  that  the  frequency  of  the  second  mode  is  matched. 

ki  =  M2K+m2)  (5) 

Equation  5  is  derived  from  | K  -  co2m|  =  0,  for  m]  =  m2  =  m,  in  which  co  is  the  angular 

frequency  to  be  matched  in  the  model.  If  the  value  of  co  is  not  readily  available,  it  can  be  easily 
determined  by  simple  testing. 

1.3.2  Mitigator 

The  representation  of  the  mitigator  is  vital  to  the  accuracy  of  the  simulation.  In  this  development, 
it  is  assumed  that  the  mitigator  is  built  of  honeycomb  and  consists  of  one  continuous  piece  of 
honeycomb  structure.  It  has  one  cylindrical  section  and  one  wedge-shaped  part  with  a  varying 
cross  section  on  the  side  (see  figure  la)  that  is  being  impacted  by  the  OBR.  A  typical  orthotropic 
crush  model  behavior  for  the  honeycomb  material  is  composed  of  three  phases  as  shown  in  figure 
3  (1 ).  These  three  phases  include  linear  elastic  loading,  volumetric  crush,  and  hardening  to  full 
compaction.  The  initial  linear  elastic  loading  phase  can  be  represented  by  an  elastic  modulus  Ee, 
yield  stress  and  strain  ay  and  ey,  respectively.  The  volumetric  crush  behavior  of  the  honeycomb 
material  model  (up  to  compaction  strain),  sc,  is  assumed  to  be  linear,  as  indicated  in  figure  3b. 

The  fluctuation  in  volumetric  strain  relations  observed  in  tests  (see  figure  3  a)  is  believed  to  have 
resulted  from  unstable  buckling  of  honeycomb  cells.  The  hardening  phase  is  assumed  to  have  a 
linear  relation  during  compaction  with  modulus  Ec.  To  model  the  mitigator,  the  relation  between 
the  contact  force  and  the  mitigator  front  deformation  and  their  progress  in  time  is  established. 
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Here,  it  is  assumed  that  as  the  wedge  section  deforms,  it  takes  the  shape  of  a  truncated  wedge,  as 
was  indicated  by  the  finite  element  analysis  (FEA)  results.  It  is  also  assumed  that  within  the 
wedge  part,  the  cross-sectional  area  varies  linearly  over  the  wedge  height.  In  the  development  of 
this  mitigator  deformation  model,  its  mass  is  herein  neglected;  its  mass  and  velocity  throughout 
the  event  are  small  compared  to  the  other  moving  parts.  The  mitigator  deformation  is  driven  in 
time  by  its  elastic  and  plastic  wave  propagation  speeds.  Since  the  elastic  propagation  speed  is 
much  higher  than  the  velocities  in  the  event,  it  is  assumed  that  elasticity-driven  events  propagate 
with  infinite  speed.  Therefore,  before  the  stress-strain  state  attains  material  yielding,  it  must  be 
homogeneous  along  the  cylinder  axis. 


AL/L0  £ 

a  b 


Figure  3.  Honeycomb  stress-strain  curve. 

Our  aim  is  to  determine  the  force  resulting  from  the  mitigator  defonnation  within  the  central 
difference  time-stepping  scheme.  Assuming  that  we  have  found  a  solution  at  time  equal  to  tn,  we 
need  to  express  the  force  at  time  tn+\  =  t„  +  At. 

When  the  OBR  hits  the  mitigator,  a  plastic  region  forms  immediately  at  the  wedge  top.  As 
indicated  in  figure  4,  the  boundary  between  this  plastic  region  and  the  rest  of  the  mitigator, 
which  is  still  elastic,  would  travel  down  the  length  of  the  mitigator  with  variable  velocity  equal 
to  the  current  plastic  wave  speed,  cpi  (2,  5).  Figure  4  also  shows  the  respective  material  states 
(density  [p],  stress  [a])  of  the  crushed  and  elastic  parts  of  the  mitigator.  The  stresses,  velocities, 
and  densities  on  both  sides  of  the  plastic  wave  front  are  further  explained  in  figure  5.  With  the 
approach  suggested  in  (2)  and  (J),  the  plastic  wave  speed  can  be  determined  as  follows. 


Plastic  wave  front 
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Figure  5.  Mitigator  deformation. 


Following  the  notation  of  figure  5,  let  the  initial  length  of  the  mitigator  be  Lq.  At  the  current 
analysis  time,  t,  the  total  deformation  of  the  mitigator  is  ALo  and  the  lengths  of  the  crushed  and 
the  elastic  parts  are  x  and  Lei,  respectively.  Since  the  yield  strain  of  the  mitigator  is  negligible 
compared  to  its  compaction  strain,  we  can  assume  that  all  mitigator  deformation,  ALo,  occurs  in 
the  compacted  region: 


Figure  6.  Mitigator  wedge  deformation  zones. 


The  mitigator  deformation  can  be  expressed  through  the  displacements  of  both  its  ends,  which 
are  equal  to  the  displacements  of  the  adjacent  masses  m2  and  m3. 

=uobr  ~umem  (7. 

The  plastic  wave  speed,  cpi,  is  equal  to  the  time  rate  of  change  in  length  of  the  crushed  portion  as 
shown  next. 
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,,  _  d(x  +  ALq)  _  1  d(AL0)  __  1  d(uOBR  uMEm)_{Vobr  V MEM ) 

c»/  —  ,  ~  '  r  ~  V0/ 

dt  ec  dt  sc  dt  sc 

Knowing  the  velocity  of  the  OBR  at  time  tn,  we  can  calculate  this  speed  and  revise  the  plastic 
wave  front  location: 

[X  front) n+\  =  {Xfr»nt)n+{Cp,)n^  (9) 

The  stresses  immediately  below  the  plastic  wave  front  (assuming  a  vertical  mitigator  with  the 
wedge  on  top,  figure  6)  are  equal  to  the  yield  stress.  Thus,  knowing  the  position  of  the  wave 
front,  one  can  easily  determine  the  area  at  the  plastic  wave  front  (Af)  and  the  corresponding 
cross-sectional  resultant  contact  force,  Fcont  as  shown  next. 


Fcont=Vyi  Af)  (10) 

The  front  area,  Af,  is  determined  from  the  plastic  wave  front  displacement.  Within  the  wedge 
portion  of  the  mitigator,  we  assume  that  there  are  two  regions:  the  plastic  zone  and  the  elastic 
zone.  As  illustrated  in  figure  6,  Lw  is  the  wedge  initial  height,  and  Ld  is  its  current  defonned 
height  to  be  determined  from  loading.  When  the  plastic  zone  length  ( Lp  )  is  less  than  the  wedge 

length  ( Lw  ),  then  the  front  area  can  be  determined  as  follows: 


Otherwise, 


At=  Ao~f~ 
Af=  Ao 


(ID 


in  which  Ao  is  the  cross-sectional  area  of  the  cylindrical  mitigator. 

Equation  10  is  assumed  to  represent  the  mitigator  force,  F,  the  contact  force  between  the  mitigator 
and  OBR,  and  the  mitigator  and  MEM  in  figure  2.  Computation  of  mitigator  force  using  equation 
10  is  believed  to  be  valid  to  a  sufficient  degree  of  accuracy  as  long  as  the  plastic  wave  front 
progresses  forward  ( cpi  >  0).  When  the  plastic  wave  motion  no  longer  propagates,  then  the  miti¬ 
gator  force  decreases  after  the  unloading  of  the  elastic  stiffness  of  the  mitigator.  This  scheme 
can  be  easily  programmed  as  a  stand-alone  routine,  as  shown  in  appendix  A. 


1.3.3  Strain  Rate  Effects 

Any  analysis  model  of  the  impact  event  would  rely  on  the  constitutive  representation  of  the 
mitigator  honeycomb.  The  associated  relationship  would  usually  be  strain  rate  dependent,  with  the 
dependency  being  very  important  for  the  accuracy  of  representation.  Significant  stress  enhance¬ 
ment  under  impact  loading  is  reported  for  most  of  the  cellular  honeycomb  materials  by  Zhao  ( 4 ). 
Since  in  most  cases  it  is  impossible  to  develop  a  closed  fonn  model  of  the  strain  rate  effect,  some 
empirical  approaches  have  been  developed  and  successfully  applied  in  this  investigation.  One 
such  example  is  a  power  law  relation  used  in  (5),  in  which  the  stress  is  expressed  as  a  function  of 
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the  quasi-static  stress-strain  relationship  and  of  the  current  strain  rate.  Experimental  data  for 
different  rates  can  be  used  to  obtain  the  parameters  of  the  power  law  relationship,  thus  expressing 
the  stress  enhancement  as  a  function  of  variable  strain  rates.  Since  no  such  experimental  data  are 
available  for  the  present  study,  the  available  quasi-static  relations  are  used,  and  a  stress  enhance¬ 
ment  coefficient  of  1.5  times  the  static  stress  is  assumed  to  model  the  crushing  stress  on  the 
honeycomb  material.  In  a  previous  finite  element  (FE)  simulation  of  the  air-gun  environment  (6), 
the  authors  have  successfully  demonstrated  the  effectiveness  of  such  a  constant  stress  enhancement 
factor  of  1.5  in  modeling  the  strain  rate  effect  on  the  honeycomb  material. 

1.4  Results  and  Discussion 


1.4.1  Validation  of  Present  Model 


To  validate  the  present  model  and  assess  its  performance,  the  predicted  responses  with  the 
present  model  are  compared  with  the  available  test  and  simulation  results.  Measured  responses 
from  the  air  gun  launch  simulation  of  the  test  item  shown  in  figure  lb,  as  well  as  the  results  from 
a  detailed  FE  model  built  to  represent  the  system  behavior  of  the  air  gun  environment,  are  used 
for  validation  and  verification  of  the  present  model. 


To  model  the  honeycomb  material,  the  static  stress-strain  curve  from  figure  3  was  used,  scaled  to 
1 .5  times  the  static  relationship  to  incorporate  the  strain  rate  effect.  This  approach  was  identical 
to  the  one  used  in  the  FE  analysis  (6).  The  total  mass  of  the  OBR  was  distributed  between  the 
two  masses  m,\  and  m2,  with  a  ratio  1:1.  To  determine  the  k\  (spring  stiffness),  the  period  of 
oscillation  (T)  of  the  OBR  was  estimated  to  be  about  1.36x10  4  s,  from  the  FE  acceleration  curve 
as  filtered  at  2500  Hz.  The  corresponding  frequency  was  determined  to  be  1/T,  or  7353  Hz,  and 
the  angular  frequency  co  equals  46,200  rad/s.  Then  from  equation  5  we  have 


<o2(m,+m2)  =  462002-3.7S  =  „ 

1  A  A 


The  values  of  all  parameters  used  in  the  present  analysis  are  now  provided: 


m  1  -  mass  of  OBR,  back:  1 .25  kg 

nil  -  mass  of  OBR,  front:  2.5  kg 

/«3  -  mass  of  MEM:  3 1 .2  kg 

k\  -  OBR  spring  stiffness:  2.0xl09  N/m 

vo  -  OBR  initial  velocity:  83.566  m/s 

AT  -  End-time  of  analysis:  0.002  s 

At-  Time-step  for  analysis:  1.0x10  6  s 

St  -  Output  time  step:  4. Ox  1 0  6  s 

sc  -  compacting  strain:  0.64 
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<jy  -  yield  stress:  55.16  MPa 

CTcrush- 1.5 -(static  yield)  =  1.5-55.16  =  82.74  MPa 

E  -  elastic  modulus:  4,060  MPa 

po  -  density  of  honeycomb:  650  kg/nr 

r  -  radius  of  mitigator:  0.05  m 

Lw  -  length  of  mitigator  wedge:  0.038  m 

Lc—  length  of  cylindrical  section:  0.219  m 


These  values,  converted  for  the  convenience  of  comparing  the  results,  were  entered  into  the  input 
file  provided  in  appendix  A.  The  analysis  was  performed  with  the  FORTRAN  (Formula  Trans¬ 
lator)  program,  and  the  corresponding  results  are  presented  and  compared  to  the  FE  analysis  and 
test  results  in  figures  7  to  16.  Figure  7  shows  the  unfiltered  acceleration  curves  of  the  OBR  top 
(attachment  point  of  mass  m2  for  present  approach).  Predicted  accelerations  filtered  at  7000  Hz, 
2500  Hz,  and  1000  Hz  are  presented  in  figures  8,  9,  and  10,  respectively.  Since  the  only  frequency 
present  in  the  current  model  is  low,  the  results  from  the  present  approach  have  been  filtered  only  in 
figure  10,  which  portrays  the  results  of  applying  a  1000-Hz  filter  to  all  three  curves.  In  figures  7 
through  9,  only  the  FE  and  test  curves  were  filtered.  As  seen  from  the  first  two  figures,  the  present 
approach  matches  the  average  of  both  the  test  and  the  FEA  curves.  Figure  9  shows  that  the  present 
model  gives  excellent  results  for  both  the  acceleration  peak  and  the  acceleration  pulse  duration. 
Figure  10,  data  filtered  at  1000  Hz,  indicates  excellent  agreement  between  the  present  and  the  FE 
approach.  Figure  1 1  shows  frequency  domain  plot  comparison  of  the  acceleration  data.  As  seen  in 
figure  1 1 ,  the  first  fundamental  frequency  of  the  present  model  matches  that  of  the  test  item.  This 
comparison  suggests  that  the  modeling  of  OBR  spring  stiffness,  k\,  based  on  test  item’s  first 
fundamental  mode,  reasonably  captures  the  dynamics  of  the  proposed  mass-spring  discrete  system 
of  the  air  gun  launch  environment. 


Time  [s]  (E-03) 

Figure  7.  Top  acceleration  of  the  OBR  unfiltered. 
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Time  [s]  (E-03) 

Figure  8.  Top  acceleration  of  the  OBR  filtered  at  7000  Hz. 


Time  [s]  (E-03) 

Figure  9.  Top  acceleration  of  the  OBR  filtered  at  2500  Hz. 
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Figure  10.  Top  acceleration  of  the  OBR  filtered  at  1000  Flz. 
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Frequency  (E+6) 

Figure  11.  Fast  Fourier  transform  (FFT)  of  the  top  acceleration  of  the  OBR. 
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OBR  velocities  for  the  test,  FE  simulation,  and  present  approach  are  compared  in  figure  12.  As 
shown  in  the  figure,  a  very  good  agreement  is  observed  between  the  FE  and  present  approach. 
Test  velocity  history  is  obtained  from  the  integration  of  test  acceleration  data.  A  similar  compari¬ 
son  of  OBR  displacement  results  for  all  three  cases  is  presented  in  figure  13.  Figure  13  shows 
that  the  present  approach  matches  the  maximum  displacement  from  the  test  very  closely.  The 
discrepancies  we  observe  to  the  right  of  the  peak  displacements  (when  the  OBR  starts  to  retract) 
can  be  attributed  to  friction  in  the  “catch  tube”  and  to  other  mechanisms  not  accounted  for  within 
the  FE  and  the  present  analysis. 


A  Test 

B  FF  Analysis 
C  Present 


Figure  14  shows  the  contact  forces  between  the  mitigator  and  the  OBR  and  between  the 
mitigator  and  the  MEM.  In  the  FE  analysis,  these  forces  differ  because  of  the  elastic  wave 
propagation  speed,  inertial  effects  in  the  mitigator,  etc.  These  phenomena  contribute 
insignificantly  to  the  dynamics  of  the  system,  especially  if  the  OBR  is  to  be  the  primary  object  of 
interest;  therefore,  they  were  neglected  within  the  present  approach.  This  simplification  makes 
the  contact  forces  at  the  two  interfaces  equal,  but  as  we  see,  they  are  very  close  to  the  average  of 
both  curves.  Furthermore,  when  we  compare  impulses  by  integrating  the  force  over  time,  we  get 
a  very  good  agreement  among  the  impulses  at  both  interfaces,  as  shown  in  figure  15. 
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Figure  13.  Top  displacement  of  the  OBR. 


ARL  GUN  SIMULATION 
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ARL  GUN  SIMULATION  -  Force  Impulse 


_A_FEA  -  OBR  Side 
B  FEA  -  MEM  Side 
C  Present  Both  Sides 


The  present  approach  provides  data  for  the  length  of  the  compacted  zone  within  the  mitigator. 
The  position  of  the  crush  front  within  the  undeformed  mitigator  (“initial”)  and  the  length  of  the 
crushed  part  (“compacted”)  are  shown  in  figure  16. 


A  Initial 
EL Compacted 
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1.5  Parameter  Variation  Sensitivity 

In  this  section,  we  illustrate  the  model’s  response  to  varying  the  parameters  of  the  system.  As 
explained  earlier,  the  representation  of  the  OBR  by  two  separate  masses  and  a  spring  aims  to 
capture  one  additional  frequency  of  the  OBR  which  is  considered  of  interest.  However,  this 
representation  requires  the  total  physical  mass  of  the  OBR  to  be  distributed  between  the  two 
representative  masses  in  the  model:  ni\  and  m2.  Furthermore,  the  connecting  spring  stiffness 
has  to  be  determined.  These  processes  introduce  additional  uncertainties  into  the  model;  this  is 
especially  true  of  the  mass  distribution  process  since  it  is  not  uniquely  defined  but  depends  on 
the  user’s  judgment.  In  general,  varying  the  m\\m2  mass  ratio  will  only  change  the  oscillation 
amplitude,  not  the  oscillation  frequency.  Figure  17  shows  the  acceleration  curves  of  masses  m\ 
and  m2  for  two  different  distribution  ratios.  First,  the  total  mass  of  the  OBR  (3.75  kg  in  this 
example)  was  distributed  to  1.25  kg  for  m\  and  2.5  kg  for  m2  (ratio  1:2).  The  corresponding 
accelerations  are  shown  with  curves  A  and  B  in  figure  17.  Curves  C  and  D  show  the  accelerations 
of  the  masses  for  a  ratio  of  2: 1  (m  1  =  2.5  kg  and  m2  =  1.25  kg).  The  frequency  of  oscillation  for  all 
four  curves  is  the  same.  Only  the  amplitude  of  oscillation  is  changed. 


Mass 

fl M1 

B  M2 
_C_M1 
D  M2 


1:2  Ratio 
1:2  Ratio 
2:1  Ratio 
2:1  Ratio 


The  stiffness  of  the  spring  connecting  the  two  OBR  masses,  k\,  is  determined  by  equation  5 
which  requires  the  determination  of  the  angular  frequency  of  the  OBR.  Since  this  would 
normally  be  detennined  experimentally,  some  error  might  be  introduced  into  the  spring  stiffness 
value.  Therefore,  it  is  interesting  to  see  how  small  variations  in  the  value  of  k\  would  affect  the 
model  output.  Effects  of  the  variation  of  k\  values  are  illustrated  in  figure  18,  which  compares 
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the  initial  acceleration  response  to  that  obtained  for  k\  values  10%  smaller  and  larger  than  the 
initial  stiffness  determined  by  equation  5.  As  seen,  the  system  behavior,  especially  in  the  area  of 
interest  (first  millisecond),  does  not  change  significantly.  Furthermore,  when  drastically 
changing  the  spring  stiffness  value,  as  shown  in  figure  19,  we  still  get  the  correct  shape  and  peak 
of  the  deceleration  pulse. 


K1  Value 

_A_2.0E+06 
JEL1  8E+06 
_C_2.2E+06 


Time  (E-03) 


Figure  18.  Model  sensitivity  to  10%  variation  of  the  value  of  kj. 


ARL  GUN  SIMULATION  M2  Acceleration  for  K1  Values 


Time  (E-03) 


Figure  19.  Model  sensitivity  to  large  variation  of  the  value  of  % 
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The  integration  time  step  is  a  very  important  parameter  of  the  central  difference  time-stepping 
scheme.  In  order  to  get  an  accurate  solution,  the  time  step  needs  to  be  smaller  than  a  certain 
value  called  the  critical  time  step.  If  this  requirement  is  not  satisfied,  the  solution  quickly 
diverges.  Therefore,  determining  the  critical  time  step  for  a  given  solution  is  very  important  to 
ensure  accuracy  in  the  results.  However,  in  most  cases  this  is  not  an  easy  task,  and  for  complex 
systems,  only  threshold  limits  of  the  critical  time  step  are  determined.  In  the  present  model, 
there  is  no  need  to  determine  the  critical  time  step  explicitly.  Since  the  solution  does  not  take  a 
significant  amount  of  run  time,  we  can  easily  re-run  it  with  a  smaller  time  step  and  compare  the 
two  solutions  to  assess  the  accuracy  of  the  results.  As  illustrated  in  figure  20,  the  initial  solution 
using  a  time  step  At  =  I  .Ox  1 0  6  s  is  identical  to  the  solution  generated  with  a  10-times-smaller 
time  step,  At  =  1.0x10  7  s.  Even  with  a  much  coarser  time  resolution  of  At  =  2.0x10  5  s,  we  get  a 
reasonably  accurate  solution.  This  shows  that  the  chosen  value  of  the  time  step  is  well  below  the 
critical  and  illustrates  the  good  overall  stability  and  robustness  of  the  solution  scheme. 


Time  Step 

_A_1  0E-06 
B  2.0E-05 
_^1  0E-07 


Finally,  to  illustrate  how  some  of  the  geometric  parameters  of  the  system  affect  the  solution 
response,  the  radius  of  the  mitigator  honeycomb  cylinder  was  varied.  The  corresponding  OBR 
acceleration  responses  are  shown  in  figure  2 1 .  As  would  be  expected,  a  smaller  radius  produced 
a  lower  acceleration  peak  and  a  wider  pulse,  whereas  a  larger  radius  produced  a  higher 
acceleration  peak  and  a  narrower  pulse. 
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ARL  GUN  SIMULATION  Mitigator  Radius  Sensitivity 
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Figure  21.  Model  sensitivity  to  10%  variation  of  the  value  of  mitigator  radius. 


1.6  Conclusion 

The  objective  of  this  effort  was  to  develop  an  overall  analytical  model  to  predict  the  transient 
response  of  a  generic  artillery  component  subjected  to  launch  simulation.  The  simulation  is 
achieved  by  impact  mitigation  techniques  in  which  the  component  is  propelled  to  a  target  inside 
a  carrier  equipped  to  measure  component  response  during  the  simulation.  A  2-degree-of- 
freedom  model  is  developed  to  simulate  the  air  gun  launch  environment  in  which  a  test  object 
mounted  on  a  projectile  is  fired  through  the  air  gun  and  decelerated  by  the  crushing  of  an 
aluminum  honeycomb  mitigator,  which  in  turn  impacts  the  MEM  before  being  stopped  at  the 
retrieving  end.  The  presented  model  achieved  a  good  prediction  of  the  period  and  peak 
acceleration  of  the  OBR. 

Sensitivity  of  parametric  study  indicates  that  the  distribution  of  masses  in  the  spring-mass  model 
of  the  test  item  does  change  the  oscillation  amplitude  without  affecting  the  system’s  frequency  of 
vibration.  A  slight  variation  of  spring  stiffness  (±10%)  does  not  change  the  shape  and  neither 
does  it  affect  the  peak  of  the  deceleration  pulse  of  the  test  item.  The  sensitivity  study  of  initial 
time  step  illustrates  the  good  overall  stability  and  robustness  of  the  proposed  solution  scheme. 
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2.  Part  II:  Air  Gun  Finite  Element  Simulation 


2.1  Review  of  Previous  Work 

Several  LS-DYNA  models  were  developed  to  simulate  an  air  gun  launch  environment.  In  this  test, 
a  test  object  is  mounted  on  a  projectile  and  fired  through  the  air  gun.  The  object  is  decelerated  by 
the  crushing  of  an  aluminum  honeycomb  mitigator  MEM  complex  before  being  stopped  at  the 
retrieving  end  ( 6 ). 

The  objective  of  the  effort  was  to  develop  an  overall  simulation  model  to  predict  the  three- 
dimensional  transient  response  of  a  generic  artillery  component  subjected  to  launch  simulation  in 
the  U.S.  Anny  Research  Laboratory’s  air  gun  facility  at  Adelphi,  Maryland.  The  simulation  is 
achieved  by  impact  mitigation  techniques  whereby  the  component  is  propelled  to  a  target  inside  a 
4-inch  carrier  equipped  to  measure  component  response  during  the  simulation.  These  models 
include  all  elements  of  the  simulation  apparatus  that  contribute  to  the  total  response  of  the 
component  under  test. 

Two  formulations  are  employed  in  the  simulations:  Lagrangian  and  arbitrary  Lagrangian/Eulerian 
(ALE).  The  aluminum  honeycomb  mitigator  undergoes  a  significant  amount  of  defonnation, 
which  could  render  a  Lagrangian  simulation  severely  unstable.  For  this  reason,  an  ALE  simulation 
is  also  considered.  The  Eulerian  method  is  more  suitable  for  problems  involved  in  severe  mesh 
distortion.  The  Lagrangian  method,  on  the  other  hand,  is  limited  in  how  much  an  element  can 
deform.  The  Lagrangian  method  is  easy  to  set  up  and  visualize  since  the  material  point  moves 
with  the  mesh.  However,  the  Eulerian  method  is  more  difficult  to  set  up  and  the  mesh  is  stationary 
so  that  material  points  are  advected  from  one  element  to  the  next. 

The  aluminum  honeycomb  mitigator  is  modeled  with  *MAT_MDIFIED_HONEYCOMB  in  the 
LS-DYNA  simulation  code.  This  is  material  model  number  126  in  the  code.  This  material  model 
is  available  for  the  Lagrangian  method  only  and  it  is  not  yet  implemented  in  the  Eulerian  solver  of 
the  LS-DYNA.  Therefore,  an  alternate  material  model  is  considered  that  can  simulate  the 
behavior  of  aluminum  honeycomb  and  is  available  in  the  LS-DYNA  Eulerian  solver.  The 
material  model  considered  for  this  case  is  *MAT_CRUSHABLE_FOAM  which  is  material  model 
number  63  in  LS-DYNA.  This  material  model  can  simulate  an  isotropic  crushable  material.  The 
*MAT_MDIFIED_HONEYCOMB  can  simulate  an  orthotropic  crushable  material.  Most  of  the 
aluminum  honeycomb  material  in  the  air  gun  simulation  is  crushed  axially.  Therefore,  the 
*MAT_CRUSHABLE_FOAM  model  is  considered  to  be  appropriate  for  this  model.  Several 
LS-DYNA  models  are  developed  to  simulate  the  air  gun  launch  environment  with  Lagrangian 
and  ALE  computational  methods.  The  Lagrangian  method  is  simpler  to  set  up,  simpler  to  post 
process,  and  require  less  computational  time.  However,  it  requires  significant  expertise  in  the 
explicit  FE  codes  to  make  the  simulation  numerically  stable.  This  is  because  of  the  significant 
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large  deformation  of  the  mitigator.  On  the  other  hand,  the  ALE  method  is  more  difficult  to  set 
up,  more  difficult  to  post  process,  and  requires  much  more  CPU  (central  processing  unit)  time. 
However,  the  ALE  simulation  is  more  suitable  for  very  large  deformation-like  flow  problems. 

Both  computational  methods  lead  to  the  same  prediction  for  the  acceleration  of  the  OBR.  An 
accurate  prediction  of  the  period  and  peak  acceleration  of  the  OBR  is  achieved  with  the  presented 
models  and  the  methods  employed. 

2.2  Finite  Element  Model  Sensitivity  Study 

It  is  very  important  that  the  FE  simulation  captures  the  contributions  of  high  frequency  oscillations 
to  the  net  acceleration  of  the  model.  This  will  allow  failure  mode  simulation  of  the  embedded 
MEMS  (micro-electrical-mechanical  systems).  Improving  the  correlation  between  test  and  FE 
results  allows  the  FE  simulation  to  become  a  tool  for  detecting/predicting  failure  modes  of 
physical  components  subject  to  an  air  gun  test.  For  this  purpose,  the  developed  FE  model  is 
further  examined  to  fully  understand  the  differences  between  the  model  and  the  experimental  data 
and  to  acquire  more  confidence  in  the  predictability  of  the  simulation.  Several  model  sensitivity 
issues  are  considered  in  detail  as  follows: 

•  double  versus  single  precision  code, 

•  time  step  size  effect, 

•  different  output  intervals, 

•  different  damping  options  in  LS-DYNA  (all  options  are  investigated), 

•  representation  of  the  OBR’s  glass  beads  with  solid  elements,  and 

•  intermittent  eigenvalue  analysis. 

Each  one  of  these  items  is  studied  in  detail,  and  the  corresponding  results  are  presented. 

2.2.1  Double  Versus  Single  Precision  Code 

The  default  LS-DYNA  code  is  a  single  precision  code.  When  the  integration  time  step  is  small 
enough  to  satisfy  the  Courant  stability  condition,  many  integration  cycles  are  possible.  Since  the 
LS-DYNA  code  is  an  explicit  code  based  on  second  order  central  difference  equations,  numerical 
error  propagation  is  possible.  The  same  model  was  run  with  single  and  double  precision  code. 

The  results  are  presented  in  figure  22.  Fast  Fourier  transfonns  (FFTs)  of  the  results  are  presented 
in  figure  23.  One  can  observe  that  both  codes  give  the  same  results. 
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Figure  22.  Acceleration  data  for  the  same  FE  model  with  two  different  precisions. 
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Figure  23.  FFT  of  the  acceleration  data. 
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2.2.2  Time  Step  Size  Effect 

Integration  time  step  size  variation  is  used  to  measure  the  stability  of  a  model.  Time  step 
variation  is  also  used  to  detect  any  anti-aliasing  in  the  collection  of  the  data  from  the  explicit  time 
integration  simulation.  The  default  integration  time  step  is  normally  taken  to  be  90%  of  the 
stability  limit  to  avoid  any  filtering  error  such  as  beating  in  the  data.  This  report  considers  two 
additional  integration  time  steps.  One  is  taken  to  be  80%  of  the  critical  time  step  and  the  other  is 
taken  to  be  40%  of  the  critical  time  step.  The  results  of  all  three  time  steps  are  depicted  in  figure 
24.  FFTs  of  the  data  are  depicted  in  figure  25.  One  can  observe  that  all  the  important  frequencies 
are  the  same  in  all  the  considered  integration  time  steps.  Therefore,  the  default  integration  time 
step  is  appropriate  for  this  problem. 
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Figure  24.  Acceleration  data  for  different  integration  time  steps. 
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Figure  25.  FFT  of  the  acceleration  data  for  different  integration  time  steps. 
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2.2.3  Different  Output  Intervals 


The  data  collection  from  the  FE  simulation  should  be  the  same  as  the  test  setup;  therefore,  the 
two  data  can  be  compared  meaningfully.  The  sampling  rate  in  the  test  is  4  x  10'6  seconds,  which 
corresponds  to  500  data  sampling  points  in  the  2-millisecond  simulation.  To  investigate  if  any 
data  are  missing  from  the  simulation  or  if  any  abnormality  exists  in  the  numerical  calculations, 
data  sampling  rates  of  10  and  20  times  the  test  sampling  rate  are  considered  (5,000  and  10,000 
points,  respectively).  In  addition,  data  are  collected  at  every  integration  time  step.  Figure  26 
depicts  the  acceleration  data  for  the  entire  sampling  rate  considered.  FFTs  of  the  data  are 
depicted  in  figure  27.  All  major  contributing  and  significant  frequencies  are  present  in  the 
lower  rate,  which  is  also  the  sampling  rate  of  the  experimental  data. 
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Figure  26.  Acceleration  data  sampled  at  different  frequencies. 
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Figure  27.  FFT  of  the  acceleration  data  sampled  at  different  frequencies. 
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2.2.4  All  Damping  Options  in  LS-DYNA  Investigated 

There  are  three  options  for  damping  in  the  LS-DYNA  code:  global,  mass,  and  stiffness  propor¬ 
tional  damping.  Global  damping  is  similar  in  fonnulation  as  the  mass  proportional  damping,  and 
therefore,  results  for  this  case  are  not  presented.  Different  damping  constants  are  considered  for 
each  damping  option.  For  the  mass  proportional  damping,  constants  of  magnitudes  100,  500, 
1000,  and  1500  are  considered.  Results  for  these  cases  are  depicted  in  figure  28.  FFTs  of  the 
acceleration  data  are  depicted  in  figure  29.  One  can  observe  that  frequencies  are  damped  signifi¬ 
cantly,  especially  the  lower  frequencies  of  the  system.  Figure  30  depicts  the  test  data  and  the 
damped  acceleration  data  from  the  simulation. 
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Figure  28.  Acceleration  data  for  simulation  with  mass  proportional  damping. 
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Figure  29.  FFT  of  the  acceleration  data  for  simulation  with  mass  proportional  damping. 
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Figure  30.  Acceleration  data  for  simulation  and  test  with  mass  proportional  damping. 


Stiffness  proportional  damping  is  considered  in  the  simulation  as  1%,  10%,  25%,  25%,  and  35% 
of  the  critical  damping.  The  acceleration  data  are  depicted  in  figure  3 1 .  FFTs  of  the  data  are 
presented  in  figure  32.  One  can  observe  that  the  low  frequencies  are  not  damped,  while  the 
higher  frequencies  are  damped.  Figure  33  depicts  accelerations  from  the  test  data,  the  simulation 
with  no  damping,  and  the  simulation  with  10%  of  the  critical  damping.  Figure  34  depicts  the 
FFT  of  these  accelerations.  One  can  observe,  again,  that  higher  order  frequencies  are  affected 
more  significantly  than  lower  order  frequencies.  In  general,  mass  proportional  damping  should 
be  used  for  low  frequency  content  and  pure  defonnation  with  no  rigid  body  motion  since  it 
damps  the  rigid  body  motion.  Stiffness  proportional  damping  is  orthogonal  to  rigid  body  motion 
and  would  not  affect  this  type  of  motion.  Stiffness  proportional  damping  is  good  for  high 
frequency  content  and  should  be  used  in  such  problems  where  attenuation  of  high  frequency 
vibration  is  desired. 
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Figure  31.  Acceleration  data  for  simulation  with  stiffness  proportional  damping. 
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Figure  33.  Acceleration  data  for  simulation  and  test  with  stiffness  proportional  damping. 
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Figure  34.  FFT  of  the  acceleration  data  for  simulation  and  test  with  stiffness  proportional  damping. 


2.2.5  The  OBR  Filled  With  Solid  Elements  to  Represent  Beads 

The  OBR  is  filled  with  beads  in  the  experiment.  However,  the  OBR  is  hollow  in  the  simulation. 
These  beads  exist  in  the  OBR  to  protect  the  accelerometer  during  tests.  There  is  no  material 
model  in  any  FE  code  to  represent  the  behavior  of  these  beads.  In  this  report,  a  model  of  the 
OBR  is  considered  that  is  filled  with  solid  elements.  Different  material  models  such  as  soil, 
foam,  viscoelastic,  and  crushable  foam  are  considered  for  the  beads  in  order  to  simulate  the 
beads’  damping  high  frequency  oscillations.  Results  of  the  simulation  are  depicted  in  figure  35. 
One  can  observe  that  high  frequency  oscillations  are  more  significant  because  of  waves  traveling 
back  and  forth  in  the  solid  element. 
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Figure  35.  Acceleration  data  for  simulation  with  filled  OBR. 
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2.2.6  Intermittent  Eigenvalues  Analysis 

To  understand  the  frequency  content  in  this  simulation,  an  intermittent  eigenvalue  analysis  is 
conducted.  The  LS-DYNA  code  allows  for  the  extraction  of  frequency  contents  during  an  impact 
simulation.  This  is  perfonned  to  relate  the  frequency  content  to  the  event  and  deformation  mode  in 
the  OBR.  The  FFT  of  the  unfiltered  acceleration  of  the  OBR  is  depicted  in  figure  36.  One  can  see 
several  important  peaks  in  the  data.  The  largest  peak  corresponds  to  the  first  fundamental  mode. 
The  next  peak  is  at  about  7400  Hz.  The  next  ones  are  at  about  8000,  14600,  21000,  23800,  27700, 
and  29500  Hz.  The  intermittent  eigenvalue  analysis  revealed  the  sources  of  these  frequencies. 

One  hundred  modes  were  extracted  from  the  analysis.  Only  the  modes  corresponding  to  the  peaks 
in  the  data  are  presented  here.  Figure  37  depicts  the  undefonned  (line  mesh)  and  the  deformed 
mode  (solid  shaded  mesh)  corresponding  to  the  peaks  in  the  frequency  spectrum  graph  (figure  36). 
It  is  clear  that  these  peaks  correspond  to  the  axial  mode  of  the  OBR.  High  frequency  content  was 
not  observed  in  the  experimental  data.  However,  high  frequency  content  is  observed  in  the 
simulation  data.  Since  the  beads  are  not  present  in  the  simulation  and  the  high  frequencies  are 
attributable  to  the  axial  deformation  modes  of  the  OBR,  one  can  conclude  that  these  high 
frequencies  in  the  simulation  results  are  not  realistic.  This  conclusion  is  based  on  the  assumption 
that  high  frequency  oscillations  in  the  axial  direction  of  the  OBR  are  absorbed  by  the  glass  beads 
surrounding  the  accelerometer. 
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Figure  36.  FFT  of  the  acceleration  data  for  simulation 
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e.  f. 

Figure  37.  Mode  shapes  for  a  fundamental  frequency. 
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2.3  Simplified  Simulation 

The  full  FE  model  runs  for  few  hours  on  a  single  processor  CPU.  For  parametric  studies,  this  run 
time  is  considered  to  be  long.  A  simplified  methodology  in  which  an  FE  model  can  run  in  few 
minutes  is  developed.  This  methodology  requires  that  a  set  of  simulations  be  run  for  each  OBR 
configuration  (see  figure  lb,  for  example)  with  the  discrete  spring-mass  model  shown  in  figure  2, 
and  extract  the  force  exerted  between  the  test  item  and  the  mitigator.  The  contact  forces  (F  in 
figure  2)  between  the  OBR  and  the  mitigator  could  subsequently  be  used  in  a  simplified  FE 
model  of  the  test  item  only.  For  a  given  contact  force  corresponding  to  a  given  OBR  weight  and 
configuration,  an  FE  mesh  of  the  OBR  (test  item)  could  be  analyzed  that  requires  a  shorter  run 
time.  Any  modification  in  the  OBR  content  that  may  be  required  because  of  the  addition  of  a 
new  electronics  package  to  be  tested  could  easily  be  modeled  in  the  simplified  FE  mesh  of  the 
OBR  for  a  given  force.  To  simulate  the  response  of  the  OBR,  the  contact  force  is  applied  as  an 
external  force  to  the  OBR  to  resist  the  momentum  of  the  OBR.  This  methodology  is  applied  to 
the  current  model  and  presented  next. 

For  the  test  item  used  in  figure  lb,  the  contact  force  between  the  mitigator  and  the  OBR  is 
presented  in  figure  38.  This  force  is  then  applied  as  an  external  force  to  the  OBR  FE  model  (the 
only  mesh  is  the  OBR)  for  a  given  initial  impact  velocity.  The  acceleration  data  are  collected 
from  the  simplified  model  and  compared  to  full  model  in  figure  39.  One  can  observe  that  they 
are  practically  the  same.  The  simplified  model  runs  in  less  than  a  minute. 
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Figure  38.  Contact  force  between  the  OBR  and  the  mitigator. 
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Figure  39.  Acceleration  data  of  the  simplified  and  full  model. 


2.4  Conclusion 

The  sensitivity  of  FE  model  parameters  affecting  the  dynamics  of  a  test  object  launched  in  an  air 
gun  test  is  discussed  in  this  part  of  the  report.  An  LS-DYNA  FE  model  previously  developed  to 
simulate  the  impact  mitigation  environment  is  used  to  investigate  the  sensitivity  of  FE  model 
parameters.  Previous  study  indicates  the  necessity  of  honeycomb  material  characterization  under 
very  high  impact  loading  for  effective  modeling  of  the  strain  rate  effects.  The  sensitivity  study 
clearly  indicates  that  the  model  is  numerically  stable  and  can  aid  in  failure  prediction  if  damping 
is  added.  A  comparison  of  FE  simulation  responses  with  various  damping  options  suggests  that 
the  mass  proportional  damping  attenuates  the  low  frequency  vibration  in  the  simulated  response 
while  the  stiffness  proportional  damping  attenuates  the  high  frequency  oscillations.  The  FE 
analysis  also  suggests  that  intermittent  eigenvalue  analyses  during  an  impact  simulation  could  be 
useful  in  identifying  the  dominant  modes  of  vibration  governing  the  dynamics  of  the  test  items 
and  verifying  the  predictive  capability  of  the  FE  model.  The  reasons  that  the  FE  model  predicted 
high  frequency  content  (above  14,000  Hz)  that  is  not  observed  in  the  test  could  include  (a)  a  lack 
of  information  for  material  characterization  of  the  honeycomb  mitigator  under  very  high  strain 
rate,  or  (b)  lack  of  a  proper  material  model  for  the  beads  which  has  significant  damping  at  high 
frequencies. 

An  efficient  procedure  using  a  simplified  decoupled  model  is  proposed  to  simulate  the  response 
of  a  test  item  launched  in  an  air  gun  test.  In  this  approach,  the  contact  force  generated  from  the 
discrete  analytical  model  is  applied  on  a  detailed  FE  mesh  of  the  test  item  for  a  given  impact 
velocity.  The  acceleration  data  obtained  by  this  procedure  compared  well  with  those  of  the  full 
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FE  analysis.  This  model  is  programmed  in  a  FORTRAN  code  that  can  produce  results  with  only 
a  few  minutes  of  run  time  which  are  comparable  to  the  FE  and  test  results.  This  procedure  can 
aid  in  the  prediction  of  the  failure  of  MEMS  in  the  OBR  with  deformable  FE  model. 
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Appendix  A.  Program  I/O  Files  and  FORTRAN  Source  Code 


Input  File 

The  input  parameters  necessary  for  the  impact  model  are  entered  in  an  ASCII  text  file  named 
“input.dat,”  which  needs  to  be  located  in  the  same  directory  as  the  program  executable.  A 
sample  input  file  with  brief  description  of  the  input  parameters  follows: 


1.25 

Ml  -  mass  of  OBR  back 

2.5 

M2  -  mass  of  OBR  front 

31.2 

M3  -  mass  of  MEM 

2.0E6 

K1  -  OBR  spring  stiffness 

83.566 

V0  -  OBR  initial  velocity 

0.002 

End  time  of  analysis 

1.0E-6 

Time  step  for  analysis 

4.0E-6 

Output  time  step 

0.64 

Compacting  strain 

82.74E3 

Yield  stress 

4 . 060E6 

Elastic  modulus 

650.0 

Density  of  honeycomb 

0.05 

Radius  of  mitigator 

0.038 

Length  of  wedge  part 

0.219 

Length  of  cylindrical  part 

Notes: 

The  input  parameters  must  be  in  consistent  units.  The  program  does  not  check  for  unit 
consistency  nor  does  it  convert  the  input  values  to  ensure  consistency.  The  units  used  in  the 
above  input  are  kilogram,  second,  meter,  Newton. 

The  program  reads  only  the  first  input  parameter  from  each  line  of  the  input  file.  Everything 
after  the  first  parameter  is  treated  as  comments.  However,  the  first  parameter  should  have  the 
appropriate  value  of  the  corresponding  variable  and  be  in  the  appropriate  numeric  format. 

The  input  parameters  should  be  entered  in  the  exact  order  as  shown  in  the  example  above.  There 
should  be  no  empty  lines  at  the  beginning  of  the  input  file. 
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Output  Files 

The  program  produces  several  output  files: 

output.dat  -  an  ASCII  file  with  all  the  output  information  in  different  columns.  There  are  12 
columns  having  the  values  of  current  analysis  time,  displacement,  velocity,  and  acceleration  of 
masses  M\,  Mo,  and  M3,  current  force  value,  and  current  plastic  front  location. 

LS_POSTdis.dat  -  an  ASCII  file  with  the  displacements  of  the  three  mass  attachment  points. 
This  file  is  in  LSTC  LS-POST  format. 

LS_POSTvel.dat  -  an  ASCII  file  with  the  velocities  of  the  three  mass  attachment  points.  This 
file  is  in  LSTC  LS-POST  format. 

LS_POSTacc.dat  -  an  ASCII  file  with  the  accelerations  of  the  three  mass  attachment  points.  This 
file  is  in  LSTC  LS-POST  format. 

LS_POSTfor.dat  -  an  ASCII  file  with  the  contact  force  between  the  mitigator  and  OBR  (MEM). 
This  file  is  in  LSTC  LS-POST  format. 
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FORTRAN  Source  Code 


The  FORTRAN  source  code  of  the  program  implementing  the  formulation  described  in  section 
“Model  Formulation”  follows  below. 


program  main 

double  precision  uml (3) , un (3) , upl (3) , udn (3) , uddn (3) , xklu, 
*  buff (10) , cpl, xfront, dxfront, Force 

c  Initialize 

xpi=3. 1415926536 

open  ( 10, file= ' input . dat ' ) 

read (10, *) xMl 

read (10, *) xM2 

read (10, *) xM3 

read (10, * ) xKl 

read (10, *) xVO 

curtime=0 . 

read (10, * ) endtime 

read ( 1 0 , * ) dt 

dt2=dt**2 

curout=0 . 

read (10, *) dtout 

mcount=0 

read (10, *) epsc 

read (10, * ) sigy 

read (10, *) Emod 

epsy=sigy/Emod 

read (10, * ) rhoO 

rhoy=rho0/ (l.-epsy) 

rhofrnt=rho0/ (l.-epsc) 

read (10, * ) radius 

read ( 1 0 , * ) xLwed 

read (10, *) xLcyl 

close (10) 

xf ront=0 . 

area0=xpi*radius**2 
Force=0 . 
c  Time  =  0 
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c  Displacement 

un (1) =0 . 
un (2 ) =0 . 
un ( 3 ) =0 . 
c  Velocity 

udn ( 1 ) =xV0 
udn (2 ) =xV0 
udn ( 3 ) =0 . 

c  Acceleration 

uddn (1) =0 . 
uddn (2 ) =0 . 
uddn ( 3 ) =0 . 

c  Displacement  at  -dt 

do  i—1 , 3 

uml (i) =-dt*udn (i) 
end  do 

c  Define  output 

open  ( 10, file= ' output . dat ' ) 
c  End  Initialization,  Define  Time  Stepping  Loop 

c 

c  Stop  when  acceleration  of  MEM  gets  zero 

c  do  while  (uddn ( 3 ) . gt . 0 . 0 ) 

c 

c  Or  stop  at  endtime 

c 

do  while  (curtime . It . endtime) 

c  Write  values  for  displacement,  velocity  and  acceleration 

if  (curtime . ge . (curout+dtout) )  then 

cur out=cur time 

write (10, 10) curtime,  (un (i)  ,  udn (i) , uddn (i) , i=l, 3) , Force, xfront 
mcount=mcount+l 
end  if 

c  Calculate  mitigator  force 

cpl= (udn (2) -udn (3) ) /epsc 
if  (cpl.gt.0.)  then 
c  Plastic  region  grows 

dxf ront=cpl*dt 
xf ront=xf ront+dxf ront 
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if  (xf ront . It . xLwed)  then 
areac=areaO*xfront/ xLwed 
else 

areac=areaO 
end  if 

Force=sigy*areac 

else 

c  Elasricity  only;  unloading 

Force=Force+ 

*  areaO*Emod* (udn (2) -udn (3) ) *dt/ (xLwed+xLcyl+un (3) -un (2) -xfront) 
end  if 

if  (Force . It . 0 . )  Force=0 . 

c  Get  new  displacement  values 

xklu=xKl* (un (1) -un (2) ) 
upl ( 1 ) =2 . *un (1 ) -uml (1) -dt2*xklu/ xMl 
upl (2 ) =2 . *un (2 ) -uml (2 ) +dt2* (xklu- Force) / xM2 
upl ( 3 ) =2 . *un ( 3 ) -uml ( 3 ) +dt2  * Force /xM3 
do  i=l , 3 

c  Velocity  and  acceleration 

udn (i)=.5* (upl (i) -uml (i) ) /dt 

uddn ( i) = (upl (i) -2.*un (i) +uml ( i ) ) /dt2 

c  Shift  displacements  in  matrix 

uml ( i ) =un ( i ) 
un ( i ) =upl ( i ) 
end  do 

c  Step  forward  in  time 

curtime=curtime+dt 
end  do 

c  Write  LS-POST  output  files 

open  (ll,file='LS  POSTdis.dat') 
write (11,20) 
do  i=l,3 
rewind (10) 

write  ( 1 1 , 30 ) i , mcount 
do  j=l, mcount 

read (10, *) cur time, (buff (k) , k=l, 10) 
write (11, 40) curtime,  buff(3*i-2) 
end  do 
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write (11, 50 ) 
end  do 
close  (11) 

open  (ll,file='LS  POSTvel.dat') 
write (11, 120) 
do  i=l , 3 
rewind (10) 

write  ( 1 1 , 130 ) i , mcount 
do  j=l, mcount 

read (10, *) cur time, (buff (k) , k=l, 10) 

write (11, 140) curtime,  buff ( 3  * i  —  1 ) 

end  do 

write (11, 50 ) 

end  do 

close  (11) 

open  ( 1 1 , f ile= ' LS_POSTacc . dat ' ) 
write (11, 220) 
do  i=l,3 
rewind (10) 

write  (11, 230) i, mcount 
do  j=l, mcount 

read (10, *) curtime, (buff (k) , k=l, 10) 

write (11, 240) curtime,  buff(3*i) 

end  do 

write (11, 50 ) 

end  do 

close  (11) 

open  (ll,file='LS  POSTfor.dat') 
write (11, 320) 
rewind (10) 

write  ( 1 1 , 330 ) mcount 
do  j=l, mcount 

read (10, *) curtime, (buff (k) , k=l, 10) 

write (11, 340) curtime,  buff (10) 

end  do 

write (11, 50 ) 

close  (11) 

close  (10) 
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10  format (11 (E12.6, ' , ' IX) ,E12.6) 

20  format  ( ' Curveplot ' , / ' ARL  GUN  SIMULATION Time ' , 

*  / ' Displacement '  ,  /  ' Mass ' ) 

30  f ormat ( ' Displ .  of  M',I1,'  #pts=',I10) 

40  format (E12 . 6, 2X, E12 . 6) 

50  format (' endcurve ' ) 

120  format  (' Curveplot ARL  GUN  SIMULATION Time ' , 

*  /' Velocity Mass ' ) 

130  f ormat (' Veloc .  of  M',I1,'  #pts=',I10) 

140  format (E12 . 6, 2X, E12 . 6) 

220  format  (' Curveplot ARL  GUN  SIMULATION Time ' , 

*  /' Acceleration Mass '  ) 

230  f ormat (' Accel .  of  M',I1,'  #pts=',I10) 

240  format (E12 . 6, 2X, E12 . 6) 

320  format  (' Curveplot ARL  GUN  SIMULATION Time ' , 

*  / ' Force '  ,  /  ' Mass ' ) 

330  format (' Contact  Force','  #pts=',I10) 

340  format (E12 . 6, 2X, E12 . 6) 
stop 
end 
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MIAMI  LAKES  FL  33014 

3  FOSTER  MILLER 

M  ROYLANCE  W ZUKAS 
195  BEAR  HILL  RD 
WALTHAM  MA  02354-1196 

1  ROM  DEVELOPMENT  CORP 
ATTN  ROMEARA 
136  SWINEBURNE  ROW 
BRICK  MARKET  PLACE 
NEWPORT  RI  02840 

1  TEXTRON  SYSTEMS 
ATTN  M  TREASURE 
1449  MIDDLESEX  ST 
LOWELL  MA  01851 

1  O  GARA  HESS  &  EISENHARDT 
ATTN  M  GILLESPIE 
9113  LESAINT  DR 
FAIRFIELD  OH  45014 

1  MILLIKEN  RSCH  CORP 

ATTN  M  MACLEOD 
PO  BOX  1926 
SPARTANBURG  SC  29303 

1  CONNEAUGHT  INDUSTRIES  INC 
ATTN  J SANTOS 
PO  BOX  1425 
COVENTRY  RI  02816 

1  ARMTEC  DEFENSE  PRODUCTS 
ATTN  S  DYER 
85  901  AVE53 
PO  BOX  848 

COACHELLA  CA  92236 


3  PACIFIC  NORTHWEST  LAB 

ATTN  M  SMITH  G  VAN  ARSDALE 
R  SHIPPELL 
PO  BOX  999 
RICHLAND  WA  99352 

5  ALLIANT  TECHSYSTEMS  INC 
ATTN  C  CANDLAND  MN1 1  2830 
C  AAKHUS  MN11  2830 
B  SEEMN11  2439 
RDOHRNMN11  2830 
D  KAMDAR  MN 1 1  2830 
5050  LINCOLN  DR 
MINNEAPOLIS  MN  55436-1097 

3  ALLIANT  TECHSYSTEMS  INC 
ATTN  N  VLAHAKUS  MN  1 1  2 1 45 
S  HAGLUND  MN  1 1  2439 
M  HISSONG  MN  1 1  2830 
4700  NATHAN  LANE  N 
PLYMOUTH  MN  55442-2512 

1  R  FIELDS 

4680  OAKCREEK  ST 
APT  206 

ORLANDO  FL  32835 

1  APPLIED  COMPOSITES 
ATTN  WGRISCH 
333  NORTH  SIXTH  ST 
ST  CHARLES  IL  60174 

1  CUSTOM  ANALYTICAL  ENG  SYS  INC 
ATTN  A  ALEXANDER 
13000  TENSOR  LANE  NE 
FLINTSTONE  MD  21530 

1  AAI  CORPORATION 

ATTN  DR  N  B  MCNELLIS 
PO  BOX  126 

HUNT  VALLEY  MD  21030-0126 

1  OFC  DEPUTY  UNDER  SEC  DEFNS 
ATTN  J  THOMPSON 
1745  JEFFERSON  DAVIS  HWY 
CRYSTAL  SQ4  STE501 
ARLINGTON  VA  22202 

3  ALLIANT  TECHSYSTEMS  INC 
ATTN  J  CONDON  E  LYNAM 
J  GERHARD 
WV01  16  STATE  RT  956 
PO  BOX  210 

ROCKET  CENTER  WV  26726-0210 
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PROJECTILE  TECHNOLOGY  INC 
515  GILES  ST 

HAVRE  DE  GRACE  MD  21078 

AEROJET  GEN  CORP 
ATTN  D  PILLASCH  T  COULTER 
C  FLYNN  D  RUBAREZUL 
M  GREINER 

1100  WEST  HOLLYVALE  ST 
AZUSA  CA  91702-0296 


2  UDLP 

ATTN  G  THOMAS  M  MACLEAN 

PO  BOX  58123 

SANTA  CLARA  CA  95052 

2  UDLP 

ATTN  R  BRYNSVOLD 
P  JANKEMS  170 
4800  EAST  RIVER  RD 
MINNEAPOLIS  MN  55421-1498 


BRIGS  COMPANY 
ATTN  J  BACKOFEN 
2668  PETERBOROUGH  ST 
HERNDON  VA  22071-2443 

ZERNOW  TECHNICAL  SERVICES 
ATTN  L  ZERNOW 
425  W  BONITA  AVE  STE  208 
SAN  DIMAS  CA  91773 

GENERAL  DYNAMICS  OTS 
FLINCHBAUGH  DIV 
ATTN  K  LINDE  T  LYNCH 
PO  BOX  127 
RED  LION  PA  17356 

GKN  WESTLAND  AEROSPACE 
ATTN  D  OLDS 
450  MURDOCK  AVE 
MERIDEN  CT  06450-8324 

SIKORSKY  AIRCRAFT 
ATTN  GJACARUSO  B  KAY 
T  CARSTENSAN 
S  GARBO  MS  S330A 
J  ADELMANN 
6900  MAIN  ST 
PO  BOX  9729 

STRATFORD  CT  06497-9729 

PRATT  &  WHITNEY 
ATTN  C  WATSON 
400  MAIN  ST  MS  114  37 
EAST  HARTFORD  CT  06108 


1  LOCKHEED  MARTIN 
SKUNK  WORKS 
ATTN  D  FORTNEY 
1011  LOCKHEED  WAY 
PALMDALE  CA  93599-2502 

1  NORTHRUP  GRUMMAN  CORP 

ELECTRONIC  SENSORS  &  SYSTEMS  DIV 
ATTN  E  SCHOCH  MS  V  16 
1745 AW  NURSERY  RD 
LINTHICUM  MD  21090 

1  GDLS  DIVISION 
ATTN  D  BARTLE 
POBOX  1901 
WARREN  MI  48090 

1  GDLS 

ATTN  MPASIK 
PO  BOX  1800 

STERLING  HEIGHTS  MI  48090-1800 

1  GDLS 

MUSKEGON  OPERATIONS 
ATTN  MSOIMAR 
76  GETTY  ST 
MUSKEGON  MI  49442 

1  GENERAL  DYNAMICS 
AMPHIBIOUS  SYS 
SURVIVABILITY  LEAD 
ATTN  G  WALKER 
991  ANNAPOLIS  WAY 
WOODBRIDGE  VA  22191 


INST  FOR  ADVANCED  TECH 
ATTN  H  FAIR  I  MCNAB 

P  SULLIVAN  S  BLESS 
W  REINECKE  C  PERSAD 
4030-2  W  BRAKER  LN 
AUSTIN  TX  78759-5329 

ATTN  M  LIN  W  WEB 
1440  N  KRAEMER  BLVD 
ANAHEIM  CA  92806 


AEROSPACE  CORP  5 

ATTN  G  HAWKINS  M4  945 
2350  E  EL  SEGUNDO  BLVD 
EL  SEGUNDO  CA  90245 

CYTEC  FIBERITE 
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ARROW  TECH  AS  SO  1 

1233  SHELBURNE  RD  STE  D8 
SOUTH  BURLINGTON  VT  05403-7700 

R  EICHELBERGER 
CONSULTANT 

409  W  CATHERINE  ST  1 

BEL  AIR  MD  21014-3613 

SAIC 

ATTN  G  CHRYSSOMALLIS 

8500  NORMANDALE  LAKE  BLVD 

SUITE  1610  1 

BLOOMINGTON  MN  55437-3828 

UCLA  MANE  DEPT  ENGR IV 
ATTN  H  T  HAHN 

LOS  ANGELES  CA  90024-1597  5 

UMASS  LOWELL 
PLASTICS  DEPT 
ATTN  N  SCHOTT 
1  UNIVERSITY  AVE 
LOWELL  MA  01854 

IIT  RESEARCH  CENTER  1 

ATTN  DROSE 

201  MILL  ST 

ROME  NY  13440-6916 

MICHIGAN  ST  UNIV 
MSM DEPT 

ATTN  RAVERILL  1 

3515  EB 

EAST  LANSING  MI  48824-1226 
PENN  STATE  UNIV 

ATTN  CBAKIS  1 

212  EARTH  ENGR 
SCIENCES  BLDG 
UNIVERSITY  PARK  PA  16802 

PENN  STATE  UNIV  3 

ATTN  RS  ENGEL 
245  HAMMOND  BLDG 
UNIVERSITY  PARK  PA  16801 

PURDUE  UNIV 

SCHOOL  OF  AERO  &  ASTRO 

ATTN  C  T  SUN  1 

W  LAFAYETTE  IN  47907-1282 


UNIV  OF  MAINE 
ADV  STR  &  COMP  LAB 
ATTN  R  LOPEZ  ANIDO 
5793  AEWC  BLDG 
ORONO  ME  04469-5793 

JOHNS  HOPKINS  UNIV 
APPLIED  PHYSICS  LAB 
ATTN  P  WIENHOLD 
11100  JOHNS  HOPKINS  RD 
LAUREL  MD  20723-6099 

UNIV  OF  DAYTON 
ATTN  J  M  WHITNEY 
COLLEGE  PARK  AVE 
DAYTON  OH  45469-0240 

UNIV  OF  DELAWARE 

CTR  FOR  COMPOSITE  MTRLS 

ATTN  J  GILLESPIE  M  SANTARE 

S  YARLAGADDA  S  ADVANI 
D  HEIDER 

201  SPENCER  LABORATORY 
NEWARK  DE  19716 

UNIV  OF  ILLINOIS  AT  URBANA  CHAMPAIGN 

DEPT  OF  MATERIALS 

SCIENCE  &  ENGINEERING 

ATTN  J  ECONOMY 

1304  WEST  GREEN  ST  115B 

URBANA  IL  61801 

MISSISSIPPI  STATE  UNIVERSITY 
DEPT  OF  AEROSPACE  ENGNRNG 
ATTN  A  J  VIZZ1NI 
MISSISSIPPI  STATE  MS  39762 

DREXEL  UNIV 
ATTN  AS  D  WANG 
32ND  &  CHESTNUT  ST 
PHILADELPHIA  PA  19104 

UNIV  OF  TEXAS  AT  AUSTIN 
CTR  FOR  ELECTROMECHANICS 
ATTN  J  PRICE  A  WALLS 
J  KITZMILLER 
10100  BURNET  RD 
AUSTIN  TX  78758-4497 

VA  POLYTECHNICAL  INST  &  STATE  UNIV 
DEPT  OF  ESM 
ATTN  MWHYER 
BLACKSBURG  VA  24061-0219 
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1  SOUTHWEST  RSCH  INST 

ENGR  &  MATL  SCIENCES  DIV 
ATTN  JRIEGEL 
6220  CULEBRA  RD 
PO  DRAWER  28510 
SAN  ANTONIO  TX  78228-0510 

1  BATELLE  NATICK  OPERATIONS 
ATTN  B  HALPIN 

313  SPEENST 
NATICK  MAO  1760 

2  US  MILITARY  ACADEMY 
ATTN  RD  HAMPTON 
BLDG  752  MAHAN  HALL 
WEST  POINT  NY  10996 

ABERDEEN  PROVING  GROUND 

1  DIRECTOR 

US  ARMY  RSCH  LABORATORY 
ATTN  AMSRD  ARL  Cl  OK  TECH  LIB 
BLDG  4600 

1  US  AMSAA 

ATTN  AMXSY  TD  P  DIETZ 
BLDG  392 

1  US  ARMY  ATC 

ATTN  CSTE  DTC  AT  AC  I  W  C  FRAZER 
BLDG  400 

1  DIRECTOR 

US  ARMY  RSCH  LABORATORY 
ATTN  AMSRD  ARL  O  AP  EG 
M  ADAMSON 
BLDG  245 

1  DIRECTOR 

US  ARMY  RSCH  LABORATORY 
ATTN  AMSRD  ARL  SL  BM  D  BELY 
BLDG  328 

3  DIRECTOR 

US  ARMY  RSCH  LABORATORY 
ATTN  AMSRD  ARL  WM  J  SMITH 

J  MCCAULEY  M  ZOLTOSKI 
BLDG  4600 

1  DIRECTOR 

US  ARMY  RSCH  LABORATORY 
ATTN  AMSRD  ARL  WM  B  (CHIEF) 
BLDG  4600 


1  DIRECTOR 

US  ARMY  RSCH  LABORATORY 
ATTN  AMSRD  ARL  WM  BA  (CHIEF) 

BLDG  4600 

2  DIRECTOR 

US  ARMY  RSCH  LABORATORY 
ATTN  AMSRD  ARL  WM  BC  P  PLOSTINS 
J  NEWILL 
BLDG  390 

7  DIRECTOR 

US  ARMY  RSCH  LABORATORY 
ATTN  AMSRD  ARL  WM  BD  P  CONROY 
BFORCH  M  LEADORE  RLIEB 
B  RICE  R  PESCE  RODRIGUEZ 
A  ZIELINSKI 
BLDG  4600 

1  DIRECTOR 

US  ARMY  RSCH  LABORATORY 
ATTN  AMSRD  ARL  WM  BD  C  LEVERITT 
BLDG  390 

1  DIRECTOR 

US  ARMY  RSCH  LABORATORY 

ATTN  AMSRD  ARL  WM  BF  S  WILKERSON 

BLDG  390 

2  DIRECTOR 

US  ARMY  RSCH  LABORATORY 
ATTN  AMSRD  ARL  WM  M  S  MCKNIGHT 
J  MCCAULEY 
BLDG  4600 

3  DIRECTOR 

US  ARMY  RSCH  LABORATORY 
ATTN  AMSRD  ARL  WM  MA  (CHIEF) 

L  GHIORSE  E  WETZEL 
BLDG  4600 

19  DIRECTOR 

US  ARMY  RSCH  LABORATORY 
ATTN  AMSRD  ARL  WM  MC  (CHIEF) 

R  BOSSOLI  ECHIN 
S  CORNELISON  D  GRANVILLE 
BHART  F  PIERCE 
ERIGAS  W  SPURGEON 
BLDG  4600 
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DIRECTOR 

US  ARMY  RSCH  LABORATORY 
ATTN  AMSRD  ARL  WM  MB  J  BENDER 

T  BOGETTI  J  BROWN  L  BURTON 
R  CARTER  K  CHO  W  DEROSSET 
G  DEWING  R  DOWDING 
W  DRYSDALE  R  EMERSON 
D  GRAY  D  HOPKINS  RKASTE 
L  KECSKES  M  MINNICINO 
B  POWERS  D  SNOHA  J  SOUTH 
M  ST  AKER  J  SWAB  JTZENG 
BLDG  4600 

DIRECTOR 

US  ARMY  RSCH  LABORATORY 
ATTN  AMSRD  ARL  WM  MD  P  DEHMER 
B  CHEESEMAN  R  DOOLEY 
G  GAZONAS  S  GHIORSE 
M  KLUSEWITZ  J  LAS AL VIA 
J  MONTGOMERY  W  ROY 
J  SANDS  DSPAGNUOLO 
S  WALSH  S  WOLF 
BLDG  4600 

DIRECTOR 

US  ARMY  RSCH  LABORATORY 
ATTN  AMSRD  ARL  WM  RP  C  SHOEMAKER 
J  BORNSTEIN 
BLDG  1121 

DIRECTOR 

US  ARMY  RSCH  LABORATORY 
ATTN  AMSRD  ARL  WM  T  B  BURNS 
BLDG  309 

DIRECTOR 

US  ARMY  RSCH  LABORATORY 
ATTN  AMSRD  ARL  WM  TA  W  GILLICH 
CHOPPEL  M  ZOLTOSKI 
BLDG  4600 

DIRECTOR 

US  ARMY  RSCH  LABORATORY 
ATTN  AMSRD  ARL  WM  TA  T  HAVEL 
JRUNYEON  M  BURKINS 
E  HORWATH  B  GOOCH 
BLDG  393 


1  DIRECTOR 

US  ARMY  RSCH  LABORATORY 
ATTN  AMSRD  ARL  WM  TC  R  COATES 
BLDG  309 

4  DIRECTOR 

US  ARMY  RSCH  LABORATORY 
ATTN  AMSRD  ARL  WM  TD  D  DANDEKAR 
M  RAFTENBERG  S  SCHOENFELD 
T  WEERASOORIYA 
BLDG  4600 

1  DIRECTOR 

US  ARMY  RSCH  LABORATORY 
ATTN  AMSRD  ARL  WM  TE  (CHIEF) 

BLDG  1116A 


DIRECTOR 

US  ARMY  RSCH  LABORATORY 
ATTN  AMSRD  ARL  WM  TB  P  BAKER 
BLDG  390 
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FOREIGN  ADDRESSES 
LTD 

R  MARTIN 
MERL 

TAMWORTH  RD 
HERTFORD  SGI 3  7DG 
UK 

CIVIL  AVIATION 
ADMIN  STRATION 
T  GOTTESMAN 
PO  BOX  8 

BEN  GURION  INTERNL  AIRPORT 

LOD  70150 

ISRAEL 

AEROSPATIALE 
S ANDRE 

A  BTE  CCRTE  MD132 
316  ROUTE  DE  BAYONNE 
TOULOUSE  31060 
FRANCE 

DRA  FORT  HALSTEAD 
P  N  JONES 

SEVEN  OAKS  KENT  TN  147BP 
UK 

SWISS  FEDERAL  ARMAMENTS 

WKS 

WLANZ 

ALLMENDSTRASSE  86 
3602  THUN 
SWITZERLAND 

DYNAMEC  RESEARCH  AB 
AKE  PERSSON 
BOX  201 

SE  151  23  SODERTALJE 
SWEDEN 

ISRAEL  INST  OF  TECHNOLOGY 
S  BODNER 

FACULTY  OF  MECHANICAL  ENGR 

HAIFA  3200 

ISRAEL 

DSTO 

WEAPONS  SYSTEMS  DIVISION 
N  BURMAN  RLLWS 
SALISBURY 

SOUTH  AUSTRALIA  5108 
AUSTRALIA 


1  DEF  RES  ESTABLISHMENT 
VALCARTIER 
A  DUPUIS 

2459  BOULEVARD  PIE  XI  NORTH 
VALCARTIER  QUEBEC 
CANADA 

PO  BOX  8800  COURCELETTE 
GOA  IRO  QUEBEC 
CANADA 

1  INSTITUT  FRANCO  ALLEMAND 
DE  RECHERCHES  DE  SAINT  LOUIS 
DE  M  GIRAUD 
RUE  DU  GENERAL 
CASSAGNOU 
BOITE  POSTALE  34 
F  68301  SAINT  LOUIS  CEDEX 
FRANCE 

1  LAB  TECHNOLOGIES  COMPOSTS 
&  POLYMERS 
J  MANSON 
DMX  LTC 

CH  105  LAUSANNE  EPFL-IMX 
SWITZERLAND 

1  TNO  DEFENSE  RESEARCH 
R IJSSELSTEIN 
ACCOUNT  DIRECTOR 
R&D  ARMEE 
PO  BOX  6006 
2600  JA  DELFT 
THE  NETHERLANDS 

1  FOA  NATL  DEFENSE  RESEARCH 
ESTAB 

DIR  DEPT  OF  WEAPONS  & 
PROTECTION 
B  JANZON 
R  HOLMLIN 
S  172  90  STOCKHOLM 
SWEDEN 

1  DEFENSE  TECH  &  PROC  AGENCY 
GROUND 
I  CREWTHER 
GENERAL  HERZOG  HAUS 
3602  THUN 
SWITZERLAND 
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1  MINISTRY  OF  DEFENCE 
RAFAEL 

ARMAMENT  DEVELOPMENT 
AUTH 

M  MAYSELESS 
PO  BOX  2250 
HAIFA  3 1021 
ISRAEL 

B  HIRSCH 

TACHKEMONY  ST  6 
NETAMUA  42611 
ISRAEL 

1  DEUTSCHE  AEROSPACE  AG 
DYNAMICS  SYSTEMS 
M  HELD 
PO  BOX  1340 

D  86523  SCHROBENHAUSEN 
GERMANY 
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